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Ground state diffusion Monte Carlo is used to investigate the binding energies and carrier probability dis¬ 
tributions of excitons, trions, and biexcitons in a variety of two-dimensional transition metal dichalcogenide 
materials. We compare these results to approximate variational calculations, as well as to analogous Monte 
Carlo calculations performed with simplified carrier interaction potentials. Our results highlight the successes 
and failures of approximate approaches as well as the physical features that determine the stability of small car¬ 
rier complexes in monolayer transition metal dichalcogenide materials. Lastly, we discuss points of agreement 
and disagreement with recent experiments. 


Atomically thin layers of crystalline transition metal 
dichalcogenides (TMDCs) have been the subject of intense in¬ 
vestigation in recent years [1,2]. As with graphene, TMDCs 
exhibit remarkable properties that originate from their quasi 
two-dimensional nature [3, 4]. However unlike graphene, 
TMDCs are direct gap semiconductors, opening up a wealth 
of potential practical applications ranging from field-effect 
transistors to photovoltaics [5, 6]. Furthermore, due to the 
lack of inversion symmetry in these single layer crystals, the 
so-called K -points on opposite corners of the two-dimensional 
hexagonal Brillouin zone are inequivalent [7-9], As a result, 
a distinct valley degree of freedom associated to states near 
these points emerges which may be manipulated and con¬ 
trolled, leading to the possibility of novel “valleytronic” ap¬ 
plications [10-13]. Lastly, carrier confinement and reduced 
dielectric screening in these materials leads to large many- 
body effects, resulting in bound state complexes of electrons 
and holes with very large binding energies [14-20]. In this 
work we focus on this latter property, providing a deeper un¬ 
derstanding of the factors that control the binding energies of 
electron-hole complexes in two-dimensional TMDCs. 

From the computational perspective, the most accurate 
means of describing excitonic properties in periodic solids 
currently available is the GW+BSE approach [21-25]. Un¬ 
fortunately, analogous fully ab initio approaches have not 
been developed for the treatment of larger electron-hole com¬ 
plexes such as trions and biexcitons [17-20]. However, sim¬ 
plified approaches have proved to be effective, building on 
well-established coarse-grained methodologies developed for 
semiconductor quantum wells and other nanostructures. An 
effective real-space electron-hole potential is combined with 
an approximate treatment of the band structure, such as an 
effective mass model or a few-band tight-binding model, to 
build the model Hamiltonian. This also has roots in the 
early discussion of the Bethe-Salpeter approach by Hanke and 
Sham, demonstrating the relationship to the phenomenologi¬ 
cal approach of Wannier [26]. 

The pioneering work of Keldysh highlighted the fact that 
screening effects in quasi two-dimensional systems are intrin¬ 
sically non-local [27]. Using a generalized Keldysh approach, 
Cudazzo et al. formulated a simple and successful theory 
for excitons in graphane [28, 29]. This approach has since 


been applied by various authors to study optical spectra as 
well as the properties of excitons, trions, and biexcitons in 
TMDCs [30-35]. These studies have produced exciton bind¬ 
ing energies and real-space structures that are in reasonable 
quantitative agreement with first principles GW+BSE calcula¬ 
tions and experiments in a variety of two-dimensional TMDC 
systems [36]. This fact is not entirely surprising for three rea¬ 
sons. First, recent ab initio calculations show that the effective 
quasiparticle interactions that emerge at the RPA level nearly 
perfectly match those used to describe the effective electron- 
hole interaction in the models mentioned above [37]. Second, 
the ab initio band structure near the K-point is well described 
by elementary two- and three-band models [12, 38]. Third, the 
spatial extent of the exciton that emerges from fully ab initio 
calculations is sufficiently large relative to the atomic scale to 
suggest that a coarse-grained Hamiltonian is justified [16]. 

Even within the simplified framework of an effective 
Hamiltonian, the exact solution of the multi-body Schrodinger 
equation for larger electron-hole complexes is challenging. 
Initial work on exciton and trion binding energies in TMDCs 
employed variational wave functions [30]. This approach has 
been used more recently and with more intricate trial wave 
functions to study biexcitons [17]. In both cases the results 
found from variational solutions of the effective few-body 
Schrodinger equation are in reasonable agreement with ex¬ 
perimental results. However, since binding energies for tri¬ 
ons and biexcitons are extracted with reference to the exciton 
binding energy, the use of variational wave functions for all 
excitonic complexes leads to binding energies that need not 
provide a lower bound to the “exact” value, and it is unclear 
how much error cancellation occurs as a result. For the trion 
binding energy, Ganchev et al. have discovered a remarkable 
exact solution, but only for the case where the full Keldysh 
effective potential is replaced with a completely logarithmic 
form that is accurate only at short range [39]. It is the pur¬ 
pose of this work to investigate the nature and accuracy of 
these approximate solutions by comparing with numerically 
exact results, and thereby to provide insights into the proper¬ 
ties of higher-order excitonic complexes in two-dimensional 
TMDCs. 

Diffusion Monte Carlo (DMC) provides a useful approach 
for studying the energetics of excitonic complexes. Briefly, 
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DMC X (eV) 

variational X 

DMC X- (meV) 

experimental X 

DMC XX (meV) 

experimental XX 

MoS 2 

0.5514 

0.54 

33.8 

43 [40], 18 [18] 

22.7 

70 [41] 

MoSe 2 

0.4778 

0.47 

28.4 

30 [42] 

17.7 


ws 2 

0.5191 

0.50 

34.0 

30 [43], 45 [44] 

23.3 

65 [43] 

WSe 2 

0.4667 

0.45 

29.5 

30 [45] 

20.0 

52 [17] 


TABLE I. Estimated exciton (X), trion (X~), and biexciton (XX) binding energies for different members of the 2D TMDC class of materials. 
Where two numbers are reported, the number on the left is the most current estimate. The statistical uncertainty in the DMC data is on the 
order of 0.1-0.3 meV. The column labeled ‘variational' refers to results based on the Keldysh form, taken from Ref. [30]. 


the DMC algorithm propagates an initial wavefunction in 
imaginary time using a Jastrow-based guiding wavefunction 
until the exact ground state wavefunction and energy is ob¬ 
tained. Technical details of our DMC calculations can be 
found in the Supplemental Material. At convergence, DMC 
yields numerically exact exciton, trion and biexciton ground- 
state energies within the confines of an effective few-body 
Schrodinger equation. Specifically, our calculations employ 
an effective mass treatment of the band structure and an 
electron-hole interaction appropriate for the two-dimensional 
TMDC family of materials, i.e. 

= 9t9jV(r u ) ( 1 ) 

i ‘ i<j 

where 

V(r) = -- [Hair/ro) - To(r/r„)]; (2) 

(ei + e 2 )h) 

in the potential above, Hq is the Struve function, To is the 
Bessel function of the second kind, and e\ and e 2 are the di¬ 
electric constants for the material above and below the TMDC 
layer; in all results presented, we use ei = e 2 = 1, relevant for 
‘ideal’ or suspended TMDC monolayers. 

In addition, DMC allows a full sampling of the square of 
the wavefunction, which can be used to extract insight into the 
structure of small bound carrier assemblies. Although DMC 
has previously been used to calculate ground-state properties 
for trions interacting with a purely logarithmic potential [39], 
to the best of our knowledge it has not been used to calculate 
trion properties with the more realistic electron-hole interac¬ 
tion above, nor has it been used to calculate the properties of 
biexcitons. It should be noted that while the present work was 
underway, a numerically exact finite temperature path integral 
Monte Carlo (PIMC) study of excitons, trions, and biexcitons 
using the full Keldysh effective potential appeared [46]. While 
we believe that DMC is a more direct method than PIMC for 
the study of what are essentially ground state properties, we 
note that the results presented here are in quantitative agree¬ 
ment with those presented earlier in Ref. [46], yielding ground 
state energies that lie below those of Ref. [46] by fractions of 
a percent. On the other hand, the goals of this work are some¬ 
what distinct from those of Ref. [46]. In particular, we focus 
on the specific physical factors that influence the delicate bal¬ 
ance of relative trion and biexciton binding energies, as well 
as the accuracy of variational approaches in light of the “ex¬ 
act” DMC results. 


In Tab. I, we report exciton, trion, and biexciton binding 
energies calculated via DMC and compare to those extracted 
in recent experiments. The DMC exciton binding energies, 
defined as E x = -E x , are only 2^1% larger than those 
obtained in previous variational calculations that employed 
a trial wavefunction of the form 'f' Tval -(r c h) ~ exp (-r ch /«). 
The radial probability distribution for the distance r c h, which 
completely determines the exciton wavefunction, is plotted in 
Fig. 1(a); we compare the variational wavefunction to results 
obtained via DMC as well as a grid-based exact diagonaliza- 
tion of the one-dimensional Schodinger equation. The sim¬ 
ple 1 s-like variational wavefunction matches the true ground- 
state wavefunction well, but does not decay rapidly enough 
for large r. As we will show, achieving a similar level of 
agreement between exact DMC and variational estimates for 
the binding energies of larger excitonic complexes is, in prin¬ 
ciple, a much more difficult task because trion and biexci- 
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FIG. 1. (a) Radial probability distributions for the distance r e h of an 
exciton in MoS 2 ; ‘ex. diag.’ refers to the result of a grid-based exact 
diagonalization. (b) Radial probability distributions for the distances 
r eh and r tt of a negative trion in MoS 2 . (c) Radial probability distri¬ 
butions for the distances r e h, r ce , and rhh of a biexciton in WSe 2 . For 
the DMC calculation, the curves for r ee and r hh coincide because the 
electron and hole effective masses are taken to be equal. 
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Trion 





Biexciton 


Keldysh 

pure In 

pure 1/r 

variational 1 

variational 2 

Keldysh 

pure In 

pure 1/r 

variational 1 

variational 2 

MoS 2 

33.8 

48.9 

1630 

26 

14 

22.7 

26.2 

2610 



MoSe 2 

28.4 

39.3 

1780 

21 

12 

17.7 

21.2 

2820 



WS 2 

34.0 

53.6 

1050 

26 

14 

23.3 

28.7 

1680 



WSe 2 

29.5 

44.9 

1120 

22 

12 

20.0 

23.9 

1780 

37 

16 


TABLE II. Comparison of trion and biexciton binding energies for several potential forms in units of meV obtained from DMC, except for the 
column labeled ‘variational’; the latter results are based on the Keldysh form and taken from Ref. [30] for trions and Ref. [17] for the WSe 2 
biexciton. Binding energies in the ‘variational 1’ column are with respect to the variational exciton binding energies, whereas those in the 
‘variational 2’ column are with respect to the exact DMC exciton binding energies. For the DMC Keldysh and pure logarithmic potentials, the 
uncertainty is of the order of 0.1-0.3 meV. For the pure Coulombic potential, the uncertainty is of the order of 10 meV. 


ton wave functions are more elaborate and their approxima¬ 
tion may in principle require many variational parameters to 
achieve a high level of accuracy. 

The DMC trion binding energies given in Tab. I are all in 
the range of 28-34 meV, which is in excellent agreement with 
current experimental estimates; however it should be noted 
that realistic substrate effects have been ignored in the present 
calculations. In Tab. II, we compare trion binding energies for 
two additional potentials: a purely logarithmic form and an 
unscreened 1/r Coulombic form. These two potentials repre¬ 
sent the asymptotic small and large r behavior, respectively, 
of (2). The purely logarithmic potential approximation has 
been employed by Ganchev et al. in their analytical treatment 
of trions in TMDCs [39] while the Coulomb potential is the 
standard form for three-dimensional semiconductors. We find 
that a purely logarithmic potential overbinds the trion and re¬ 
sults in binding energies about 50% larger than those reported 
by experiments. Unsurprisingly, the pure Coulombic potential 
vastly overbinds the complex, resulting in binding energies 
that are 30-50 times too large and with a different ordering 
than is the case for the full potential (2), which is material 
dependent. Coulombic binding energies would of course be 
reduced with the inclusion of a static dielectric constant sig¬ 
nificantly larger than unity. 

Our trion binding energies are about 30% larger than those 
calculated variationally. The two-parameter variational trial 
wavefunction used in the trion calculations was [30] 

v I / T,varOe 1 h, r e2 h) ~ exp (-r e , h /a - r e2 h/b) + [a <-> b}, (3) 

a form inspired by Chandrasekhar’s treatment of the hydro¬ 
gen anion [47]. Although Fig. 1(b) shows that this optimized 
variational wavefunction reproduces p(r e h) almost exactly, the 
variational form does not capture the electron-electron repul¬ 
sion properly because it lacks any explicit r c - r e , correlation 
terms. The peak of the electron-electron distribution is at too 
small a radius, which results in over-estimating the electron- 
electron repulsion and underestimating the trion binding en¬ 
ergy, as seen in Tab. II. Nonetheless, the level of agreement is 
surprisingly good given the simplicity of the variational wave 
function employed in Ref. [30]. Furthermore, by treating the 
exciton and trion on equal footing with physically similar vari¬ 
ational wavefunctions, a fortuitous cancellation of total energy 
errors leads to binding energies which are quite close to the 
exact results (‘variational 1’ column in Tab. II); referencing 


the variational trion energy to the exact DMC exciton energy 
(‘variational 2’ column in Tab. II) leads to a significant under¬ 
estimation of the binding energy, albeit one that is a genuine 
lower bound. 

Finally, in Tabs. I and II we report biexciton binding ener¬ 
gies E xx = 2 E x - E xx , and in Fig. 1(c) we compare car¬ 
rier probability distributions obtained from DMC and from a 
recent six-parameter variational calculation [17]. Other than 
the /'hh distance, which is accurately predicted, the variational 
wavefunction is a bit too compact, leading to a total variational 
energy which is slightly too high. When referenced to the ex¬ 
act DMC exciton energy, we note that the variational biexci¬ 
ton binding energy is actually quite accurate - only 3.5 meV 
(about 15%) too small. When referenced to the less accu¬ 
rate variational exciton energy, the mis-matched cancellation 
of errors is such that the biexciton binding energy is slightly 
overestimated. 
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FIG. 2. (a) Radial probability distributions for the distances r e h and 
r ec in the trion (X~) and biexciton (XX) using a Keldysh form for 
the inter-carrier potential, (c) The same r ee distributions as in panel 
(a), but the relocated repulsive weight is shaded. (b),(d) The same 
as in panels (a),(c) but using a Coulombic form for the inter-carrier 
potential. 
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In our DMC calculations, we find that the full potential 
with dielectric screening (2) yields binding energies that are 
smaller than either of the other potentials considered, and that 
are smaller than experimental estimates by 60-70%. Most 
significantly, exact DMC calculations show that the binding 
energies for biexcitons are significantly smaller than that of 
trions. This fact, which has been noted in recent PIMC cal¬ 
culations as well [46], disagrees with recent experimental es¬ 
timates and is at odds with expectations that emerge from the 
standard case of pure Coulombic interactions. Whereas in the 
latter 1/r case biexcitons are more strongly bound than trions 
by a factor of about 1.6, in the purely logarithmic case the 
situation is reversed and trions are more strongly bound than 
biexcitons. Interestingly, we find that for realistically parame¬ 
terized Keldysh potentials, the biexciton binding energies are 
slightly smaller than those found with the purely logarithmic 
potential, despite the latter being a presumably “weaker” po¬ 
tential; this highlights the subtle balance of energies involved 
in the formation of the biexciton. We note in passing that 
the binding energies for biexcitons obtained with the logarith¬ 
mic potential are significantly closer to the full Keldysh re¬ 
sults than they are for trions, suggesting that the short-range 
approximation of Ganchev el al. may be even better for biex¬ 
citons. This result is consistent with the smaller real-space 
structure of the biexciton seen by comparing Figs. 1(b) and 
(c). 

To gain deeper insight into this balance of energies, we con¬ 
sider the electron-hole and electron-electron distributions for 
trions and biexcitons obtained with the Keldysh and Coulom¬ 
bic potentials, plotted in Fig. 2(a), (b). Suppose that for a given 
potential, the attractive electron-hole probability profiles were 
identical for the trion and the biexciton, and the repulsive 
electron-electron (hole-hole) profiles were also identical for 
the trion and the biexciton. Then elementary arguments us¬ 
ing the definition of the trion and biexciton binding energies, 
along with the pairwise additive potential, show that the biex¬ 
citon binding energy would be exactly twice the trion binding 
energy, E xx /E x - 2. Any deviations from this ideal ratio are 
due to relative differences in the attractive and repulsive prob¬ 
ability distributions as the second hole is added to the negative 
trion. 

Instead, biexciton-to-trion binding energy ratios of less than 
2 are observed for both the screened interaction (2) and the 
Coulomb interaction. In both cases, about 9% of the total 
weight in the p(r e h) (attractive) profile is relocated from long r 
to short r. More significantly, a much larger fraction of the to¬ 
tal weight in the p{r ee ) (repulsive) is relocated to short r, lead¬ 
ing to a reduced biexciton-to-trion ratio. Specifically, for the 
screened potential (2), about 31 % of the weight is relocated, 
which leads to a biexciton binding energy that is smaller than 


the trion binding energy; for the Coulomb potential, only 26% 
of the weight is relocated, and the biexciton binding energy 
remains larger than the trion binding energy. The relocated 
repulsive area is shaded for both potentials in Fig. 2(c), (d). 
From an energetic standpoint, this more notable change in the 
repulsive profile occurs because in the trion, there is no reason 
for the like charges to be physically close in space. However, 
in the biexciton, the complex can achieve stabilizing electron- 
hole interactions by having the like charges closer together in 
space. 

A final question that may be raised concerns the qualita¬ 
tive difference between biexcitonic stability as found by DMC 
calculations and that extracted from experiments. As men¬ 
tioned above, experimentally reported biexciton binding en¬ 
ergies significantly exceed experimentally determined trion 
binding energies, and are about a factor of two or more larger 
than calculated DMC values. Since our DMC values are exact 
within the confines of the effective mass and effective poten¬ 
tial models, one possibility is that these model ingredients are 
oversimplified, and need to be amended. Although we have 
neglected screening from the substrate and surrounding en¬ 
vironment, the results of Ref. [46] suggest that the biexciton 
binding energy will remain less than the trion binding energy, 
in so far as these latter screening effects can be captured in 
the effective Keldysh potential. In particular, it is unclear if 
the assumption of effective pair-wise additive interactions is a 
good one for larger excitonic complexes; perhaps three-body 
or higher-order interactions are needed. On the other hand, ex¬ 
perimental determination of biexciton binding energies in the 
TMDCs is quite difficult and involves both assumptions of the 
nature of spectral signals as well as extrapolations. Clearly 
future work should be devoted to addressing this interesting 
discrepancy between theory and experiment. 

Note added- Since this work was completed, a preprint 
has appeared which uses high-accuracy stochastic variational 
Monte Carlo to calculate the properties of excitons, trions, 
and biexcitons in monolayer TMDCs [48]. The results are in 
agreement with the present manuscript and the authors further 
speculate as to the origin of the biexciton discrepancy noted 
above. 
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Supplemental Material: 

Binding energies and spatial structures of small carrier complexes 
in monolayer transition metal dichalcogenides via diffusion Monte Carlo 


1 


COMPUTATIONAL DETAILS 


In this section we outline the computational approach and model used to investigate excitonic complexes. Our technique of 
choice is diffusion Monte Carlo (DMC). Because DMC is such a widely used approach, we do not give a detailed account of the 
method, and refer the reader to more complete technical discussions [1-3]. 

DMC calculations use the imaginary time Schrodinger equation along with a guiding wavefunction 'P T to project out excited 
states from an initial wavefunction <l>(f = 0), propagating it in imaginary time until the true ground state wavefunction \f/Q is 
found. If we define the mixed probability /(R, t) = 'F X (R)®(R, t) [4, 5], where R contains the spatial coordinates of all quantum 
particles, then the equation of motion for /(R, t) derived from the imaginary time Schrodinger equation is 

- 0 + Vr • [v(R)/(R,0] + (£l(R, 0 - £ref)/(R, t), (SI) 

where 

v(R, 0 = ^ 1 (R)VrT t (R), (S2) 

£l(R) = v Ft(R)~ 1 //'1 , t(R), (S3) 

and E re f is an arbitrary energy offset. A solution to the importance-sampled imaginary-time Schrodinger equation is then sampled 
using approximate Greens functions that result in the drift-diffusion motion and the branching action [1]. For the systems we 
consider, the exact ground state wave function is nodeless, so DMC yields exact ground state energies and a sampling of the 
exact ground state wavefunctions. 

The guiding wavefunction used in this work is of the form 'Ft(R) = e y(R) , which contains the Jastrow factor introduced 
in Ref. [6] adapted specifically to the potential (2). The Jastrow term contains two-body electron-hole and electron-electron 
(hole-hole) terms 

Uch(r) = c\r 2 ln(r)e~ C2r - c 3 r (l - e~ C2r2 ), (S4) 


M ee (r) = c 4 r 2 ln(r)e C5, ~. 


(S5) 


The constants c i = m e mh/2(m e + mu) and c 4 = -m e /4 (m e ,h are the effective masses of the carriers) were chosen to satisfy the 
logarithmic analogue of the Kato cusp conditions; the other constants were optimized via unreweighted variance minimization 
to improve the efficiency of the Monte Carlo sampling [7, 8]. A blocking method is used to gauge the correlation timescales 
for the energy estimates, and yields accurate standard deviations for the final average [9]. Energy estimates were obtained for 
calculations with At e {0.01,0.03,0.1], and then extrapolated to zero timestep. All reported DMC probability distributions were 
sampled from forward-walking DMC calculations [10] with the optimal guiding wave function described by Eqs. (S4)-(S5), and 
At = 0.01. A forward walking time of 300 a.u. was used for calculations employing the Keldysh potential; that time was 30 
a.u. for calculations using the Coulomb potential. 

Screening lengths and effective masses for all materials studied were determined in Ref. [11]. Potentials of the type (2) were 
first discussed by Keldysh [12] and have been used to treat excitonic properties in quasi two-dimensional materials in several 
recent studies [11, 13-19]. The potential (2) behaves as 1/r at large r/r 0 , but diverges more weakly as ln(r/r (l ) near r — 0. The 
crossover point is related to the screening distance r (l = 2nxm, where xm is the two-dimensional polarizability of the TMDC 
layer. For computational efficiency and consistency with past variational calculations, we use a modified form of the effective 
potential (2), given by 


V'{r) = - 

h) 


In 


r + r 0 


+ (y - In 2)e 


-r/r 0 


(S6) 


where y is Euler’s constant [13]. Calculations with the unaltered potential (2) typically result in exciton ground-state energies 
that are merely 2-3 meV lower than those produced with (2). 


[1] R. J. Needs, M. D. Towler, N. D. Drummond, and R L. Rios, J. Phys. Condens. Matter 22, 023201 (2010). 







2 


[2] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 

[3] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001). 

[4] R. C. Grimm and R. G. Storer, J. Comp. Phys. 7, 134 (1971). 

[5] M. H. Kalos, D. Levesque, and L. Verlet, Phys. Rev. A 9, 2178 (1974). 

[6] B. Ganchev, N. Drummond, I. Aliener, and V. Fal’ko, Phys. Rev. Lett. 114, 107401 (2015). 

[7] C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60, 1719 (1988). 

[8] N. D. Drummond and R. J. Needs, Phys. Rev. B 72, 085124 (2005). 

[9] H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. 91, 461 (1989). 

[10] R. N. Barnett, P. J. Reynolds, and W. A. Lester, J. Comp. Phys. 96, 258 (1991). 

[11] T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Phys. Rev. B 88, 045318 (2013). 

[12] L. V. Keldysh, JETP Lett. 29, 658 (1979). 

[13] P. Cudazzo, I. V. Tokatly, and A. Rubio, Phys. Rev. B 84, 085406 (2011). 

[14] P. Cudazzo, C. Attaccalite, I. V. Tokatly, and A. Rubio, Phys. Rev. Lett. 104, 226804 (2010). 

[15] G. Berghauser and E. Malic, Phys. Rev. B 89, 125309 (2014). 

[16] C. J. Zhang, H. N. Wang, W. M. Chan, C. Manolatou, and F. Rana, Phys. Rev. B 89, 205436 (2014). 

[17] M. M. Fogler, L. V. Butov, and K. S. Novoselov, Nat. Commun. 5, 4555 (2014). 

[18] L. Wang, A. Kutana, and B. I. Yakobson, Annalen Der Physik 526, L7 (2014). 

[19] F. Wu, F. Qu, and A. H. MacDonald, Phys. Rev. B 91, 075310 (2015). 


